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ABSTRACT 

Potential-density pair basis sets can be used for highly efficient TV-body simulation 
codes, but they suffer from a lack of versatility, i.e. a basis set has to be constructed for 
each different class of stellar system. We present numerical techniques for generating a 
biorthonormal potential-density pair basis set which has a general specified pair as its 
lowest order member. We go on to demonstrate how the set can be used to construct 
TV-body equilibria, which we then evolve using an TV-body code which calculates forces 
using the basis set. 



Key words: methods: numerical 
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1 INTRODUCTION 

JV-body simulations of collisionless stellar systems generally 
require the evaluation of a smooth gravitational potential 
from a particle distribution containing fewer particles than 
stars in the real system by many orders of magnitude. One 
method for doing this originally proposed by Clutton-Brock 
(1972,1973) involves the use of a set of basis functions which 
constitute potential-density pairs. The method has subse- 
quently been used for stability studies of galactic equilibria 
(e.g. Allen, Palmer & Papaloizou 1990, Hernquist & Ostriker 
1992, Earn & Sellwood 1995). The method expands the po- 
tential in terms of the specified basis set. The expansion 
coefficients are evaluated as integrals from the particle dis- 
tribution which can be regarded as providing a Monte-Carlo 
estimate for them. The method relies upon the ability to ac- 
curately represent the initial equilibrium distribution and 
its subsequent evolution with the major contribution to the 
potential coming from the first few members of the basis set. 
These vary on the global length scale of the entire system 
and therefore have their expansion coefficients determined 
most accurately. 

Provided the specified basis set retains these advantages 
during the evolution of the system, the particles used to 
sample the distribution will evolve in a smooth potential, 
with the added computational advantage that the A-body 
problem is reduced to N 1-body problems together with the 
evaluation of a relatively small number of expansion coeffi- 
cients. 

Earn and Sellwood (1995) performed a stability study of 
a thin disc in which they made a comparison of normal mode 
growth rates obtained using simulation techniques with val- 
ues obtained from a linear perturbation analysis of Kalnajs 



(1976). They concluded that the basis function method is 
the optimal one for the study of the linear growth of insta- 
bilities in collisionless equilibria and therefore we suppose 
small departures from the equilibrium in general. However, 
they make the comment that the method is very specific. In 
particular the basis set has to be tailored to the particular 
density distribution under investigation, which is usually a 
difficult task and this may create problems when the density 
distribution changes significantly in a simulation. 

Up to now A-body simulations have been performed us- 
ing basis sets which can be specified analytically and which 
are based on spherical or thin disc systems, though basis sets 
for more general axisymmetric systems have been developed 
(see Robijn & Earn 1996, Earn 1996). In this paper we in- 
vestigate the possibility of numerically constructing a basis 
set which can in principle be tailored to any density dis- 
tribution, although here we focus on axisymmetric systems. 
The idea is to construct the natural biorthonormal basis for 
the distribution such that the lowest order member corre- 
sponds to a specified system. The basis set is constructed by 
finding the eigenfunctions of a Fredholm integral equation 
numerically. These are used to provide the potentials, densi- 
ties and forces of the basis set on a grid. As an illustration of 
the method we have determined the basis sets appropriate 
to a class of perfect oblate spheroids and Kuzmin-Toomre 
discs. The basis set generation code requires a few hours 
CPU time to generate 40 basis potential-density pairs and 
the associated force functions on a dedicated workstation. 

As we are interested in general distributions, the equi- 
libria cannot in general be populated with particles by sam- 
pling a distribution function as this will be unknown. Instead 
we use an orbit based method (see Schwarzschild 1979). This 
populates a time-average distribution of orbits which sam- 
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pies the mass distribution using a quadratic optimisation 
technique. The subsequent departures from the specified 
equilibrium and relaxation of the iV-body system obtained 
with evolution under the action of forces computed using the 
basis set are studied. In general the equilibrium virial ratio 
is maintained to within two percent using 20K particles and 
30 basis functions. Hernquist and Ostriker (1992) noted that 
discrete particle noise causes particle energies to undergo a 
random walk and thus leads to a relaxation of the system 
away from its original configuration. This effect is also stud- 
ied and found to lead to a relaxation time of roughly order 
N crossing times. Thus we have found that the basis sets we 
have constructed can be used effectively in simulations of the 
axisymmetric equilibria we have considered. They provide a 
promising tool for problems of stability and other problems 
where the deviation from the initial configuration is not too 
large such as finding the response to weak external forcing 
and we shall report on these in the near future. 

In section 2 of this paper we describe the basic theory 
behind the construction of a customised basis set. In section 
3 we present details of the numerical construction of basis 
sets for two contrasting distributions by two different meth- 
ods. In section 4 we go on to show how the basis sets can 
be employed to construct iV-body representations of equi- 
librium stellar systems. In section 5 we evolve the equilibria 
using an TV-body simulation code which uses the appropri- 
ate basis set for gravitational force calculation and perform 
various tests. Finally in section 6 we discuss our results, and 
outline future work. 



2 GENERATING THE BASIS SET 
2.1 Theory 

Any potential-density pair (<I>o,po) satisfies Poisson's equa- 
tion, 



V 2 $ = inpo, 



(1) 



where we have adopted units such that the gravitational 
constant G = 1. The solution is given by the Poisson integral 



(2) 



This and subsequent similar integrals are taken over the vol- 
ume occupied by the mass distribution. 
Equation (^) can be rewritten as: 



$o(r)= K(r,r')w(r')$ (r')dV' 



where the weight w(r) is given by 

/ \ Po(r) 
™ (r) = -i^(r)' 

and the kernel K(r, r') by 



(3) 



(4) 



(•») 



Note that for an equilibrium (<3?o,Po) pair, w(r) is positive 
definite which is guaranteed for a density that is positive 
everywhere. This is necessary for the method given below to 
work. From (^|) it follows that 



Co(r) = / *r(r,r%(r')d^', 



where 



£o(r) = ^/iu(r)$ (r), 
and the symmetric kernel 
K(r, r') = \J w{v)w{v')K{r, r'). 



(6) 



(7) 



(8) 



This means that £o is an eigenfunction of the eigenvalue 
problem defined by the Fredholm integral equation (for a 
detailed explanation of the properties of this type of eigen- 
value problem that we use below see Courant & Hilbert, 
1955, ch. 3), 



(9) 



with eigenvalue Ao = 1. Given the eigenfunctions {£k} we 
can calculate the via (^), and the corresponding {pk} 

are given by 



AfePfe(r) = -$ fc (r)w(r). 



(10) 



The variational properties of the eigenvalue problem, taken 
together with the fact that K (r, r') is always positive, guar- 
antee that the lowest order eigenfunction which is associ- 
ated with the largest eigenvalue does not change sign and is 
therefore £o- The largest eigenvalue is Ao = 1 and Afc — > 
for k —> oo. We define an inner product 



$ k (r) Pl (r)dV. 



(11) 



Then, when suitably normalised, the eigenfunctions obey a 
biorthonormality relation 



($k,pi) = S M . 



(12) 



The completeness properties of eigenfunction expansions of 
the type we consider described in Courant & Hilbert (1955) 
chapters 3 and 5, indicate that the set is able to represent 
uniformly and absolutely any potential $ derived from the 
Poisson integral applied to a continuous density distribution 
that vanishes sufficiently rapidly at large distances. Thus 



$(r) = y^c fc $fc(r), 



where 



c fe = <$ fc ,p) = ($,Pfc>. 



(13) 



(14) 



which is all that is needed for application to iV-body calcu- 
lations. Also we have the formal expansion for the density 



P(r) = ^c fc p fc (r), 



(15) 



Thus given a density p we can obtain the corresponding $, 
and vice versa. 



2.2 Numerical solution of the eigenvalue problem 

Generally equation ([]) will not have an analytic solution 
and therefore numerical methods are required. We work in 
cylindrical polar coordinates (R, z, <f>). Although in principle 
there are no symmetry restrictions on the equilibrium pair 
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($0, po), we shall in this paper assume them to be indepen- 
dent of <f). Then the equilibrium is axisymmetric and the 
basis set can be separated into subsets with harmonic az- 
imuthal dependence associated with a particular azimuthal 
mode number m. Thus we seek a representation of a given 
potential-density pair as follows; 

p(R, z, <j>) = pim(R, z){c lm cos ra0 + d lm sinm^}, 



E 



$lm(R, z){c lm COS 7770 + d lm sin 7770}. (16) 



That is, for a given m we require basis pairs ($; m , pi m )e 



I = 0, 1, 2, . . . which satisfy (|S|), with a different weight func- 
tion w m — —pom/$>om for each value of m. 

We wish the lowest order pair ($00, Poo) to correspond 
to the axisymmetric equilibrium. Thus we take poo = Po to 
be the equilibrium density. 

To derive suitable {pom : m > 1} we begin by not- 
ing that Poisson's equation for the axisymmetric equilibrium 
distribution is 

9 2 $oo . 1 9$oo . d 2 $ 



dR 2 



1 d$oo 
R dR 



+ 



dz 



= 47rp 



(17) 



Differentiating with respect to 7? gives 



d 2 
OR 2 



/9£oo\ 
V dR J 



1 d 



RdR\ dR 



d 2 /ggoo 

dz 2 \ 
d$> 00 



R 2 dR 



= 47T 



dR 
dpoo 



dR 



(18) 



but this implies 



\dR 2 



_ 1A 
+ RdR + 



+ 



1 d 2 

R 2 d<j> 2 



)( 



di> 



dR 



4tt 



00 i<b 

e 



dp, 



00 icj> 



dR 



(19) 



Hence ( a ^oo , ^fjf 1 ) is a potential-density pair with the 
4> dependence required of a dipole (m = 1) term in 
(^). It corresponds to a small translation of the system. 
Appropriate higher order ($om,pom) can be obtained af- 
ter further differentiation with respect to R. One finds 
that a potential-density pair generated from ($00, Poo) is 
(R m (^$00) ,.R™ (^rpoo)) e m t where y = R 2 . Func- 
tions derived in this way can be used to specify a pair 
($Om,pom) provided these are both of uniform sign as is 
the case for the examples considered below. But note that 
the procedure can also be applied to every member of a basis 
set for 777 = to generate a basis set for any non-zero m. 

We present two alternative methods for solving the 
eigenvalue problem given by equation (^). The first, method 
a), involves discretising the eigenvalue problem onto a grid 
and calculating the integral of (^) via direct summation, 
whilst the second, method b) uses a basis set built from 
standard orthogonal polynomials to represent the eigen- 
value problem, the integration being calculated via Gaus- 
sian quadrature. Method a) has the advantage of in general 
being more accurate for a given amount of computer time, 
but method b) has more flexibility. 



2.3 Basis set generation via direct summation 

We seek to numerically evaluate (^) using the {pome" 71 *} as 
the density source terms. Because of the simple (j> depen- 



dence the <j> integration can be expressed as a function of 
elliptic integrals, leaving 

J Rt=0 J z , = . x ^/(R + R>) 2 + (z-z>) 2 

x p 0m {R',z')R'dR'dz', (20) 



where 



4RR' 



(R + R') 2 + (z - z') 2 



and 

F m {k 2 ) =4 



r/2 e 2im<t> c 



(1-fc 2 COS 2 (f)) 1 / 2 ' 



(21) 



(22) 



The (R 1 , z') integration is done by finite summation approx- 
imation on an 77 x 77 (Ri,Zj) grid, i.e. 



where the summation convention is assumed, and 



(23) 



[$0m]ij = ®0m{Ri,Zj) 



y/(Ri + R P ) 2 + 

4RiR p 



[Son 



dRrj 



dz a 



{Ri + R p ) 2 + { Zj - Zq y 

p0m(Rp, z q )RpdR p dz q , 

1 



R v 



(24) 



The contribution to the integral of a cell in which r = r' is 
approximated by setting the value of the kernel in the cell 
to be the average value of the kernel in the four nearest 
neighbouring cells. 

In practice, especially when equilibria are not bounded 
in space, it is necessary to transform coordinates from the 
(R, z) to the (u, v) plane, say, where the exact choice of 
transformation is chosen to achieve maximum accuracy. This 
leads to an extra Jacobian factor in [Som]„„, but does not 
affect the overall technique. 

Once this has been done the matrices are manipulated 
in the same manner as the original functions were in the 
theoretical case, giving the discrete analogue of (0), 



Aom [^On 



K ijpq [£om]p9 ■ 



(25) 



By mapping the 77 x 77 grid to a vector with n 2 elements, 
equation ( |25| ) becomes a symmetric matrix eigenvalue prob- 
lem which can be be solved, for example, using a QL factori- 
sation algorithm. Thus the values of the first 71 x 77 eigenfunc- 
tions on the coordinate grid are obtained, which can then 
be used for interpolation. 

The CPU time required to generate 40 basis functions 
and the associated forces is a few hours on a Sun Sparc20 
workstation, which is small in comparison to the time re- 
quired to conduct large long-term simulations. This suggests 
that it may be possible to generate a new basis set mid- 
way through the simulation, if the system under study has 
evolved away from its initial state to such an extent that the 
old basis set is no longer able to accurately represent it. 
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2.4 Basis set generation via Legendre polynomials 
and Gaussian quadrature 

The eigenvalue problem denned by can be expressed in 
matrix form using a polynomial basis. To do this we suppose 
a coordinate transformation that maps the cylindrical coor- 
dinates (R, z) to (u, v). This is such that the infinite domain 
< R < oo, — oo < z < oo is mapped into the finite domain 
ui < u < U2, Vi < v < D2. For examples see below. 
We now introduce the orthonormal polynomials 

Vw (l*u,fiv) = N-, 1/2 Pi (fiu) Pi> (jiv) , (26) 

for I = 0,1,2,...,/' = 0,1,2,..., where P denotes the Leg- 
endre polynomial, with both fi u — ^ 2 "~ 2 1 1 ~ l " 2 ' ) and fi v = 

^(v^-vi^ both lying in (—1, 1). Here the normalization 
constant N w = (v 2 - Vi)(w 2 -u 1 )/((2l+ l)(2l' + 1)). These 
polynomials form an orthonormal basis over the (u, v) do- 
main. 

To obtain a matrix eigenvalue problem we first expand 
(e.g. for m = 0) the eigenfunction £ko in the form 



, w ' (27) 

W V I' 

where |J| gives the Jacobian such that 

RdRdz -> \ J\dudv. (28) 
Then the inner product 

(£ko,£k'o) = -2ff y^Ofc,H>Ofc/,ij/ = -27ra fe -a fe /, (29) 

w 

may be written as a scalar product of the column vectors a^ 
and a fc / whose components are the expansion coefficients of 
£fc and £ fc / . 

Substituting the expansion (E7|) into the eigenvalue 
problem defined through (^) and using the orthonormality 
of the polynomial basis expansion, we obtain the exactly 
equivalent matrix eigenvalue problem 

A fc a fc =M-a fc , (30) 
where the symmetric matrix M has elements given by 



My 



K{r,r')Vi{u,v)Vj{u',v')y/\jy\J'\ 
X d(j>dudvdu dv' . 



(31) 



Here we have adopted a mapping of the pairs of subscripts 
associated with the basis functions into single integers i and 
j and |J'| = |J(u',»')|. 

Up until now the matrix representation of the eigen- 
value problem is exact. Approximations to it are obtained by 
truncating the polynomial basis such that it becomes finite. 
Thus we adopt I = 0, 1, 2..., L max -1,1' = 0, 1, 2..., L max - 1 
giving L max functions in each of the coordinate directions 
or I/ max in total. Solution of the now finite matrix eigen- 
value problem, which can be accomplished with the QL al- 
gorithm as above, requires evaluation of the matrix elements 
given by (J3TJ) . To do this requires evaluation of multiple in- 
tegrals. We performed these using Gaussian quadrature in 
each of the (u, v) coordinate directions using L max weights. 
This requires evaluation of the polynomial basis set at the 
coordinate points corresponding to the zeros of -PL max (Atu) 



and fL max (/iu). These points also provide a natural grid on 
which to evaluate the eigenfunctions once their expansion 
coefficients have been found from the solution of the eigen- 
value problem. For practical evaluation, the kernel K(r, r') 
was of necessity softened for r = r' by replacing \z — z'\ by 
0.05i? rather than zero. Tests showed results to be indepen- 
dent of the details of the softening. 

We found that a successful solution to the eigenvalue 
problem and construction of the eigenfunctions for the ex- 
amples given below could be obtained along these lines using 
black box software found in Press et al. (1996). Compara- 
ble results to those obtained with the grid based method, 
for similar computational overhead, for the low order func- 
tions and performance in iV-body codes was obtained for 
i max = 30. Forces could be obtained as described below 
or by simple low order interpolation and numerical differ- 
entiation of the potentials. To conduct simulations using 
a basis set it is obviously necessary to possess the forces 
associated with the basis potentials. If the potentials were 
known analytically, the simplest way to calculate the forces 
would be via direct differentiation. But because the poten- 
tials are known only on a coordinate grid or via polynomial 
approximation it is more accurate to obtain the forces via 
integration of the densities over the force kernel, i.e. for the 
density function pi m (R, z)e m ^, 

[.' _ r ) e im W '-<t>) 



F im (r) 



r 1 !!'<:> 



-pim(R', z')dV', 



(32) 



using the same numerical integration algorithm that was 
used to calculate the potentials. Hence we obtain the set of 
forces {F; m } which correspond to the potential-density pair 
basis set {$ im ,p im }. 



3 TWO EXAMPLE BASIS SETS 

As illustrations of the technique basis sets have been con- 
structed for two contrasting distributions, the perfect oblate 
spheroid (POS), and a Kuzmin-Toomre (KT) disc with 
a Gaussian vertical density profile. Henceforth all results 
shown were achieved using basis sets generated by method 
a), though similar results were obtained using functions gen- 
erated via method b). 



3.1 Perfect oblate spheroid 

The perfect oblate spheroidal potential (POSP) and its cor- 
responding density are best described in some form of pro- 
late spheroidal coordinates. Details and properties of this 
distribution can be found in de Zeeuw (1985), and an ana- 
lytic but non orthonormal set of functions was constructed 
for application to it by Syer (1995). The choice of coordi- 
nates made here for numerical integration purposes was 

1 



■ tan titans, 



R z = 



1 



(tan u — e )(e 



S [arctan(e), 7r/2), 

G [— arctan(e), arctan(e) 



(33) 



where (it, v) are orthogonal coordinates in the (R, z) plane. 
The parameter e describes the flatness of the distribution 
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Figure 1. Some example (^i m ,pi ri 
(0,0),(1,0),(5,0),(0, 1). 



) of the perfect oblate spheroid with e = 5. From top to bottom they are (l,m) 



(e.g. e — gives a spherical distribution, and the system 
flattens with increasing e). Lines of constant it give ellipses, 
lines of constant v give hyperbolae orthogonal to the con- 
stant u ellipses and we are only concerned with a finite do- 
main in the (it, v) plane. 

The perfect oblate spheroidal density and potential are 
given by: 



p(u,v) = 



$(u, v) 



e(l + e 2 

e(it tan it — v tan v) 
tan 2 u — tan 2 v 



Poo, 



(34) 



(35) 



Because the equilibrium potential is known analytically, the 
numerically obtained {$om} an d their corresponding forces 
can be directly compared with the exact values and the accu- 
racy of the numerical basis function calculation determined. 
We find that method a) using a 50 x 50 (u, v) grid is sufficient 
to produce results with a maximum relative error in the po- 
tential of less than 0.2% and in the force of 0.7% over 90% 
of the mass of the system, whilst method b) using 30 x 30 
Gaussian weights has roughly double the error of the direct 
method. This level of accuracy is more than sufficient for the 



purpose of iV-body simulations where particle discreteness 
of the iV-body system is the dominant source of error. 

Four of the potential-density pair basis functions for an 
e = 5 perfect oblate spheroid are shown in Figure [lj. From 
top to bottom they are the original distribution, the first 
radially oscillatory function, the first z oscillatory function, 
and the first dipole (m = 1) function. Having constructed 
a basis set, we should ensure that the functions are capable 
of representing the kind of perturbations away from equi- 
librium that we wish to investigate, using only the first few 
members of the set. We can illustrate this by representing 
the following potential (following Syer f995), 



$ pert (r) =$ 1 (r, ei )+a$ 2 (r,e 2 ) 



(36) 



where $1 is the POSP that the basis set was based on, with 
eccentricity ei, $2 is the POSP of a different eccentricity 
e 2 , and a is some positive number less than unity. (The in- 
tegrations involved in evaluating the expansion coefficients 
were calculated numerically using the same coordinate grid 
on which the functions were originally generated) . The con- 
vergence of the series is demonstrated graphically in Figure 
|i| for the values ei = 5, e 2 = 2, and a = 0.3. The left 
plots show successive approximations of the potential and 
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Figure 2. a) shows the basis representation of the potential given in equation ( |36[ ) in the z = plane using upto Q ma x basis functions 
and b) shows the associated relative error, whilst c) and d) are the same plots, for the radial force expansion (for z = 0). 



radial force in the z — plane using one, four, and eight 
basis functions and the right plots show the corresponding 
relative error (where we have defined the relative error of a 
quantity f, say, to be / apP rox//exact - 1). 

The high level of accuracy in the expansion representa- 
tion of the perturbed distribution gives us confidence in the 
utility of our basis set. (The error in the force expansion at 
small R is not a problem in practice as the magnitude of the 
radial force tends to zero at small R and it affects such a 
small fraction of the total mass of the system). 

3.2 Kuzmin-Toomre disc 

The Kuzmin-Toomre disc profile with a gaussian vertical 
density profile (e.g. Sellwood 1994) is given by 



(2tt) 1 /2 0o 



2*§ 



Poo, 



with 



M 



2na 2 V + a 2 



-3/2 



(37) 



(38) 



Again for numerical purposes we make a coordinate trans- 
formation in which the domain of interest is finite, 

R = atan(u), 
z = v2zq tan(u), 
u e [0,tt/2), 



v e (-7r/2,7r/2). (39) 

The potential of the density given by (|37]) is not known an- 
alytically. However, using a finite difference method on a 
equally spaced 50 x 50 (u, v) grid, the potential can be cal- 
culated to a specified accuracy. This can be compared with 
the potential obtained via the approximation ( p3| ) calculated 
on the same (u, v) grid. The maximum relative error is less 
than 1% over 80% of the mass of the disc using method a), 
which is acceptable, particularly if the particle discs we use 
in simulations are truncated at this mass level. Similar ac- 
curacy was obtained with method b) with L ma x = 30. Some 
members of the set with zo = 0.1, a = 1, M = 1 are shown 
in Figure [| From top to bottom they are the original dis- 
tribution, the first radially oscillatory function, the first z 
oscillatory function, and the first dipole (m = 1) function. 
We note that this system constitutes a severe test of the 
technique, because of the extreme differences in behaviour 
in the R and z directions. 



4 CONSTRUCTION OF EQUILIBRIUM 
MODELS 

The task is to assign positions and velocities to N discrete 
particles so that they represent a steady state system de- 
scribed by the potential-density pair ($oo,Poo)- 
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P 
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z 



Figure 3. Some example (<t>i m ,P; m ) of the KT disc with zq 
(0,0), (1,0), (5,0), (0,1). 

4.1 Orbit representation 

The approach we follow is to write the density distribution 
of the particles as a weighted sum of contributions from M 
different orbits, that is 

M 

p(r,t) =^2m k 6(r - r k (t)), (40) 
fc=i 

where each v k (t) describes a different orbit labelled with the 
index k. The total number of particles in each orbit is m& if 
each particle has unit mass (as is assumed throughout this 
paper). 



4.2 Fitting the prescribed density 

Using the first Q ma x axisymmetric members of the basis set 
constructed earlier, for any density p(r, t) we can construct 
the quantity 




Z 



0.1, a = X,M = 1. From top to bottom they are (l,m) = 

Qmax-l Qmax-1 

D{t) = (*oo- ^2 c ioM*io,Poo - ^2 c 'o(t)pw}, (41) 
where 

M 

Cio(i) = ($io(r),p(r,t)) = -^m fc $ !0 (r fc (t)). (42) 

k = l 

Using the biorthonormality relation (|l^), D can be written 
as: 

Qmax-l 

D = l-2c o+ ^ c «' ( 43 ) 
!=0 

We see that D is minimised when coo = 1, and c;o = 
for I > 1. 

Because D is non-negative and is zero only when 
p(r, t) = poo(r, i), we can regard it, in an appropriate norm, 
as a measure of the 'distance' that p(r, i) is 'away from' 
poo(r). 
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Further, as we are seeking a steady state solution which 
should have a time independent density, we can take a long- 
time average to remove the explicit t dependence in the co- 
efficients ci (t). Then we find 



D = ($00 - c io $ io ,poo- X! Cl °P l °)' 



(44) 



where each coefficient is replaced by its time average denoted 
with an overline 



cm = lim • 

t— >oo 



/ * cio(t)dt 



t 



In practice we have found that stable time averages can be 
obtained after following an orbit for a time of order a hun- 
dred crossing times for the cases we have considered. 

The task is to assign a different mass to each orbit to 
make the distribution of the particles match poo as closely as 
possible, in a time averaged sense. This can be regarded as 
being equivalent to minimising D subject to the constraint 
that mfc > 0. Using the time averaged coefficients, D can be 
written as: 



-2 

c i0 , 



D = 1 - 2c 00 + Yl 

1=0 

with 

M 

cm = - y^m fc $i (r fc ), 

k=l 

and 

M 

cf = ^2 m fc mp*io(r fc )$jo(r p ). 



(45) 



(46) 



(47) 



We may write this in the equivalent form 

D = 1 + 2m • ¥oo + m A m, (48) 

where the vector components [m] k = mk, and [*oo] k = 
3>oo (fk), with the matrix element 



minimises the deviation from the desired density distribu- 
tion whilst simultaneously using surplus degrees of freedom 
manifest in the degeneracy of the pure problem to minimise 
the profit function. 

We note that Tremaine, Henon & Lynden-Bell (1986) 
argue that in the presence of small scale mixing, all func- 
tions of the form — J C(F)dxdv, where F is the distribution 
function and C is any convex function, increase during colli- 
sionless relaxation, though in general this conjecture is false 
(Sridhar 1987). This suggests that we seek a profit function 
of the form X)fc=i C( m *0- In this context, if phase space 
is divided up into cells of unit volume, and the above in- 
tegral approximated by a sum, the {m^} may be regarded 
as the discrete equivalent of the distribution function. The 
quadratic nature of the original D motivates us to construct 
a profit function which is likewise quadratic as the nature of 
the optimisation problem then remains the same. Hence we 
have taken the profit function to be of the form 



El rn± 
\ Pk 



(50) 



where the {pk} are specified statistical weights. In practice 
a reasonable choice for these is such as to give each orbit 
equal weight, e.g. 



Pk = M. 



(51) 



Proceeding as outlined above, we form a new function D' to 
be minimised such that 



D' — D + fj,S, 



(52) 



where p is a small positive constant whose magnitude re- 
flects the degree of smoothness imposed on the solution. 
By solving this quadratic optimisation problem the required 
{mk} for an equilibrium solution can be calculated. We use 
the two phase quadratic programming method of Gill et al. 
(1986) which is implemented in the NAG library routine 
E04NCF. On each orbit k we then place n particles, where 
n is proportional to m*,. By placing the particles at equal 
time intervals along each orbit we ensure a smooth distribu- 
tion in orbital phase. 



A fe p = ^2 $;o(rfc)$;o(r p ), k,p = 1, 



,M. 



(49) 



Thus, to calculate m in such a way as to provide the closest 
fit to the specified equilibrium, we have a quadratic optimi- 
sation problem, namely to minimise D subject to the con- 
straint m k > 0, k = 1, M. 

However there is a further complication, as the number 
of eigenfunctions used is (usually considerably) smaller than 
the number of orbits, making the problem ill-conditioned, 
and the particular method we use to minimise D gives a 
solution where very few occupation numbers are non-zero. 
This ill-conditioning can be removed by the use of a profit 
function. By adding a function to D which is small when the 
distribution of orbits is smooth, the degeneracy of the solu- 
tion is removed, and a wider spread of occupation numbers 
is achieved (For examples, see Richstone & Tremaine 1988, 
Merritt 1993, Syer & Tremaine 1996). We can interpret the 
solution to the modified minimisation problem as one which 



4.3 Application to the perfect oblate spheroid 

The POSP, being a Stackel potential, is completely inte- 
grable. Each orbit has three conserved actions associated 
with it making the construction of an M-orbit family on 
which to base our equilibrium easier. But the method can 
be applied to more general mass distributions, as will be 
seen in the next section, provided an M-orbit family can be 
constructed which adequately samples phase space. 

As a consequence of integrability, a general orbit in the 
POSP separates into librations between fixed values of the 
u and v coordinates, so that orbits can be classified by the 
triple (u+,U-,v) where u+,U- are the extrema in the u co- 
ordinate and ±u the extrema in the v coordinate. In terms 
of a thin disc model we can think of oscillations in u and v as 
corresponding to radial and vertical oscillations respectively. 
Here we seek an equilibrium distribution which consists of 
stars on infinitesimally thin short-axis tube orbits (a "shell 
orbit" model). These orbits have no libration in the u direc- 
tion. Thus, in the thin disc limit they correspond to circular 
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Figure 4. Comparison of the analytic density (solid line) with 
that of an ./V-body sample (dashed line) for an e = 1 perfect 
oblate spheroid. 

orbits with superposed small amplitude vertical oscillations. 
Bishop (1987) tackled this problem via a numerical solution 
of the corresponding density/distribution function integral 
equation. 

To set up an orbit family we divide the (it, v) plane into 
a grid which evenly samples the mass of the system and 
then allocate an orbit to each grid cell which has its turning 
points randomly allocated within the cell. These orbits are 
then integrated in the POSP for such a time that the <E>;o 
are determined to a sufficient accuracy. 

Using an orbit family of 1000 orbits together with 40 
of the <3>;o and applying no smoothing (i.e. p, = 0), a good 
JV-body realisation of an equilibrium e = 1 perfect oblate 
spheroid (in the sense of a small value of D') is found 
but only around 100 orbits are occupied. However, setting 
p = 10~ 4 increases the number of occupied orbits to around 
800 whilst having little effect on the value of D' . A contour 
plot comparing the analytic density with the basis set repre- 
sentation of the iV-body distribution (using 100K particles) 
is shown in Figure ^ (where the contours are logarithmically 
spaced) . 

4.4 Application to the KT disc 

The orbits we use to build our KT disc from are taken by 
sampling one of the distribution functions used in Sellwood 
& Merritt (1994). They used an analytic distribution func- 
tion for the planar motion (derived in Kalnajs 1976). They 
assigned z velocities by first integrating the one-dimensional 
Jeans equation for the vertical motion 

18^) = 
p az 

(where F z is the vertical component of the gravitational ac- 
celeration) to find the velocity dispersion a m (z), and then by 
assuming the vertical velocity distribution to be a Gaussian 
of width a w (z). We generated our orbit family by sampling 
the distribution function in two stages. First we sampled val- 
ues of 7? and z from the density distribution p(R,z). Then 
for each particle position we sampled the velocity component 



Figure 5. Comparison of the analytic density (solid line) with 
that of an Af-body sample (dashed line) for zq = 0.1 KT disc. 

perpendicular to the vertical z axis using the Kalnajs distri- 
bution function. Finally the vertical velocity component v z 
was sampled from the Gaussian of width a w (z), the latter 
having been determined from (^). Positions and velocities 
determined in the above manner were used to provide the 
initial conditions for the orbit family. 

We adopted the planar distribution function which is 
independent of angular momentum, and hence isotropic 
(tuk = 3 in the notation of Sellwood & Merritt 1994), and 
take the z scale height zo = 0.1 As before an orbit fam- 
ily of 1000 orbits and 40 of the <£>;o were used, and setting 
p = 10 -3 resulted in around 900 of these orbits being oc- 
cupied. A contour plot comparing the analytic density with 
the basis set representation of the iV-body distribution (us- 
ing 100K particles) is shown in Figure |E| (again with loga- 
rithmically spaced contours). 



5 EQUILIBRIUM SIMULATIONS 

Our TV-body simulation code uses a standard second-order 
time-centred leapfrog integration scheme. The forces on each 
particle are calculated from the basis set expansion 

F(r fc ) = -^c, m F lm (r fc ), (54) 

l ,m 

where the fourier coefficients c; m are calculated from the set 
of particle positions {r^} via 

JV 

Giro = -}^lm(r k ), (55) 

To assess the quality of our simulations, we conduct 
energy and angular momentum conservation tests and cal- 
culate the relaxation rate for the POS and KT disc equilib- 
ria generated in the previous section. We use 20K particles, 
30 of the basis functions, and a step-size of approximately 
0.004 times the dynamical time tdyn of each system, where 
we have defined td yn to be the time for a particle at the 
half-mass radius of the system to complete a quarter of one 
circular orbit around the system. We adopt the appropriate 
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tdyn as the unit of time in both cases. The results presented 
here were produced using functions generated via method 
a), though method b) gives similar results. 

To suppress any instabilities we impose axisymmetry 
and reflection symmetry in the z — plane, by omitting 
all m > 0, and z antisymmetric functions respectively from 
our basis set expansions. We also ran the systems for several 
time units before beginning the tests to allow any initial fluc- 
tuations to phase mix away. This ensures that we are testing 
the quality of the code, rather than discreteness effects in 
the initial equilibria. We give results for the simulations over 
a period of ten dynamical times, though we ran the systems 
for many more dynamical times to ensure that the results 
given here are typical. 

Sellwood & Merritt (1994) advocated the use of a quiet 
start in their iV-body stability studies of KT discs. One fea- 
ture of their quiet start is that they allocated particles in 
groups obeying certain symmetries (e.g. in the z direction 
and in azimuthal segments) . This ensures that the net linear 
and angular momentum of the system is zero and that low 
azimuthal mode numbers can be made to be absent from 
the initial density distribution. It also ensures that any sub- 
sequent rapid growth of such anm/0 component can be 
ascribed to instabilities. However, since here we are testing 
axisymmetric equilibria simulations and effectively impose 
such symmetries on the system ab initio by excluding all 
z antisymmetric and non-axisymmetric functions from the 
basis expansions such a quiet start becomes redundant (ex- 
cept in the total momentum conservation tests). In fact we 
find that in this situation using a quiet start may raise the 
noise level and relaxation rate, since, for example, allocating 
particles in pairs, at phase positions (x, y, ±z, v x ,v y , ±u z ) ef- 
fectively halves the number of particles the expansion can 
'see' and thus raises the noise level by a factor of \/2 since 
noise scales as 1/ \/~N . 



5.1 Virial ratio 

The scalar virial theorem (see, e.g. Binney & Tremaine 1987) 
tells us that 



-I = 2KE total + PE total , 
where 



JV 

i=5> 

1=1 

-P -Etotai 



2 



(56) 



(57) 



5 EM 

JV 

2 ^ Ixi - 

-It,' 



(58) 



In exact steady state equilibrium with an infinite num- 
ber of basis functions, the quantity I should be identically 
zero. Thus the behaviour of — 2KE to t&i/PE t otai is an indi- 
cation of how close our jV-body systems are to true equilib- 
rium. As shown in Figures [] and Q the virial ratio in both 




4 6 
number of dynamical times 



Figure 6. Evolution of virial ratio of POS equilibrium simulation. 




0.990 
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number of dynamical times 



Figure 7. Evolution of virial ratio of KT disc equilibrium simu- 
lation. 



systems is maintained at unity ±0.02 which is consistent 
with the statistical error we would expect with 20K parti- 
cles. 



5.2 Total energy conservation 

The total energy of the system is given by 

-Etotai = KE totSl i + P-Etotai) (59) 

The fluctuation in Etotai is less than 0.1 % in our simulations. 

5.3 Relaxation 

In a coUisionless system, individual particle energies should 
be conserved. In our simulations, however, due to discrete- 
ness noise individual particle energies will tend to random 
walk slightly with time - i.e. the system relaxes away from 
the initial coUisionless equilibrium. This effect is shown in 
Figures ^ and [] where we plot the relative change in energy 
of the TV particles over the duration of the test of lOtdyn 
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Figure 8. Individual particle energy deviations for the POS equi- 
librium simulation. 



Figure 10. Relaxation of POS equilibrium simulation. 
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Figure 9. Individual particle energy deviations for the KT disc 
equilibrium simulation. 



against their initial energies. The vertical band-like struc- 
ture observed in these scatter plots is found in similar sim- 
ulations described in Vine & Sigurdsson (1998). They state 
that there will be regions where the discrepancy between the 
value of the potential given by the basis expansion and the 
true value of the potential is greater than average (because 
of the expansion truncation), and that the bands contain 
the particles which occupy these regions. 
To quantify the relaxation we define 



.4 



1\ 

-y 



E h (t)-E h (t ) 



E k {t ) 



(60) 



where Ei(t) is the energy of particle i at time t. We note 
that if the energies of individual particles undergo a random 
walk, A as well as each individual term in the summation is 
expected to increase linearly with time. 

In Figures [lo] and |ll] we can see that indeed A grows 
linearly in time as expected. However, the rate of increase is 
different for the two systems. There are at least two possible 
explanations for this. The first is that the definition of tdyn 
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Figure 11. Relaxation of KT equilibrium simulation. 

is fairly arbitrary and is only an indicator of the dynamical 
time-scale of the two systems, hence we may not be exactly 
comparing like with like. The second, which may explain 
why the relaxation rate of the POS simulation is smaller 
than that of the KT disc is due to the different nature of 
the orbits in the two systems. The KT disc model contains 
many low angular momentum, radial orbits which pass close 
to the centre of the system. It is likely that such orbits will 
suffer more relaxation on average as the centre of the system 
is subject to proportionately more noise disruption. With 
the POS "shell orbit" model fewer particles pass near the 
centre of the system, and the orbits are more isolated from 
each other in phase space, so we might expect the relaxation 
rate in this model to be smaller. 

For the POS simulation the relaxation time is of order 
10 5 times tdyn, whilst for the KT disc it is of the order 10 4 
times tdyn, hence we can run simulations for hundreds of 
dynamical times before becoming too concerned about the 
effects of relaxation. The relaxation rate for runs using the 
basis set generated by method b) with L ma x = 30 has about 
a 10% larger relaxation rate, because of the lower accuracy 
of the method b) functions. 
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Runs with different numbers of particles indicate that 
the relaxation rate is inversely proportional to the number 
of particles used, as we would expect. 

We also briefly study the dependence of the relaxation 
rate on the number of basis functions used. A simulation us- 
ing only the zeroth order basis function has a much reduced 
relaxation rate (by around a factor of 10), but a simulation 
using the first five z and axisymmetric functions has much 
the same relaxation rate as the original simulation which 
employs thirty. This effect can be understood from Wein- 
berg (1993), which predicts that the relaxation rate of an 
./V-body simulation code in which small scale particle fluc- 
tuations have been suppressed will be dominated by large 
scale modes which are generated by particle noise and ampli- 
fied by self-gravity. In our results the runs with five functions 
adequately represent these dominant modes and additional 
functions contribute little to the rate of relaxation. 



5.4 Linear and angular momentum 

Any simulation code in which the forces are not calculated 
via pair-wise summation will violate momentum conserva- 
tion to some degree. The drift in linear momentum can be 
compensated for by recentreing the coordinate grid on the 
centre of mass of the system. By imposing suitable symme- 
tries on the initial particle distribution its total angular mo- 
mentum L = (L x , L y ,L z ) is set to zero. During the simula- 
tions L z is conserved exactly for individual particles (which 
is imposed by the exclusion of non-axisymmetric forces), 
hence total L z is also conserved exactly. Whilst individual 
L x and L y should not be conserved, the total L x and L y are 
seen to be conserved exactly in our simulations. 



6 CONCLUSIONS AND FURTHER WORK 

We have shown that is possible to numerically construct 
a potential-density pair basis set which is customised to a 
general equilibrium distribution, and demonstrated how the 
set can represent deviations from the original distribution 
with a high level of accuracy. 

We have demonstrated how the set can be used to con- 
struct an TV-body equilibrium sample of a density distribu- 
tion without explicit knowledge of the integrals of motion of 
the system and without explicit reference to a known dis- 
tribution function using an orbit based method. We have 
developed an iV-body simulation code which calculates par- 
ticle forces using the basis set, and have demonstrated the 
accuracy of the code with regards to conservation of integrals 
of motion and relaxation effects, and described the effect of 
changing the number of particles and basis functions used. 

Because the basis construction is completely general, 
we can apply the technique to systems for which no suitable 
analytic basis set is known. Since the CPU time required to 
generate a new basis set is not prohibitive, it should also 
prove possible to follow the evolution of a given system even 
if it has evolved significantly away from its original configu- 
ration by generating a new basis set from the current density 
distribution. 

Because the computational cost of the simulation code 
scales linearly with N, large numbers of particles can be 



employed which will reduce the effect of particle noise re- 
sulting in more accurate simulations. There has been much 
discussion on the relative merits of different JV-body simu- 
lation methods (e.g. Hernquist & Barnes 1990, Hernquist & 
Ostriker 1992, Earn & Sellwood 1995, Sellwood 1997) and 
it seems clear that for certain classes of problem, namely 
detailed studies that involve not too large deviations from 
an initial distribution, a basis set expansion TV-body code is 
one of the most competitive. 

We comment further that the method can be applied 
to problems involving two (or more) distinct mass distribu- 
tions, for example a galactic disc embedded in a dark matter 
halo system. This involves combining the kind of axisym- 
metric basis set described in this paper to represent the disc 
distribution with a similarly generated spherical basis set 
to represent the halo distribution. In addition to stability 
studies, problems involving the warping of a galactic disc 
embedded in a halo, or the effects of an infalling satellite on 
the system can be studied. We hope to report on such work 
in the near future. 
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